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Abstract. This paper presents a simple and elementary proof of Butcher’s 
theorem on the order conditions of Runge-Kutta methods. It is based on 
a recursive definition of rooted trees and avoids combinatorial tools such as 
labelings and Faa di Bruno’s formula. This strictly recursive approach can 
easily and elegantly be implemented using modern computer algebra systems 
like Mathematica for automatically generating the order conditions. The full, 
but short source code is presented and applied to some instructive examples. 


1. Introduction 

A first step towards the construction of Runge-Kutta methods is the calculation 
of the order conditions that the coefficients have to obey. In the old days they 
were obtained by expanding the error term in a Taylor series by hand, a procedure 
which for higher orders sooner or later runs into difficulties because of the largely 
increasing combinatorial complexity. It was a major break-through when Butcher 

published 1963 his result of systematically describing order conditions by rooted 
trees. The proof of this result has evolved very much in meantime, mainly under the 
influence of Butcher’s later work j^] and the contributions of Hairer and Wanner 

1^]. In this paper we will present a simple and elementary proof of Butcher’s 
theorem by using very consequently the recursive structure of rooted trees. This 
way we avoid lengthy calculations of combinatorial coefficients, the use of tree- 
labelings, or Faa di Bruno’s formula as in [^, ^. Our proof is very similar in spirit 
to the presentation of B-series by Hairer in Chapter 2 of his lecture notes . 

As early as 1976 Jenks posed the problem of automatically generating order 
conditions for Runge-Kutta methods using computer algebra systems, but no replies 
were received. In 1988 Keiper of Wolfram Research concluded that the method of 
automatically calculating Taylor expansions by brute force was bound to be very 
inefficient. Naturally he turned to the elegant results of Butcher’s. Utilizing them, 
he wrote the Mathematica package Butcher.m, which has been available as part 
of the standard distribution of Mathematica since then. This package was later 
considerably improved by Sofroniou and offers a lot of sophisticated tools. 

While teaching the simple proof of Butcher’s result in a first course on numerical 
ODEs, the author realized that the underlying recursive structure could also be 
exploited for a simple and elegant computer implementation. This approach differs 
from the work of Sofroniou in various respects. We will present the full source code 
in Mathematica and some applications. 
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2. Runge-Kutta methods 

The Runge-Kutta methods are one-step discretizations of initial-value problems 
for systems of d ordinary differential equations, 

x'=f{t,x), x{to)=xo, 

where the right-hand side / : [to, T] x C K x ^ is assumed to be sufficiently 
smooth. The continuous evolution x{t) = of the initial-value problem is 

approximated in steps of length r by ~ This discrete evolution 

\^t+T,tx ig defined as an approximation of the integral-equation representation 

rt+T 

J /(cr, dcr 

by appropriate quadrature formulas: 

\ 

t -f CiT, X + t'^ Oijkj . 

1=1 J 

The vectors ki G i = 1,..., s, are called stages, s is the stage number. Following 
the standard notation, we collect the coefficients of the method into a matrix and 
two vectors 

A = (a,,),, G R*^^ b = (5i,..., bsf G R^ c = (ci,..., Csf G R^ 

The method is explicit, if A is strictly lower triangular. The method has order 
p G N, if the error term expands to 

In terms of the Taylor expansion of the error <I> — 4', the vanishing of all lower order 
terms in r just dehnes the conditions which have to be satisfied by the coefficients 
A, b and c of a Runge-Kutta method. 

If we choose Ci = i shown |^, that there is no loss of generality 

in considering autonomous systems only, i.e., those with no dependence of / on t. 
Doing so, the expressions and are likewise independent of t. We 

will write and '^'^x for short, calling them the flow and the discrete flow of the 
continuous resp. the discrete system. 


( 1 ) 


= X + biki, ki = f 


i=l 


3. Elementary differentials and rooted trees 


The Taylor expansions of both the phase flow ^'^x and the discrete flow '^'^x are 
linear combinations of elementary differentials like 


nffj'fj) 


dxidxjdxk dxi dxm 

.0 h i. m 


We will use the short multilinear notation of the left hand side for the rest of the 
paper. 

An elementary differential can be expressed uniquely by the structure of how 
the subterms enter the multilinear maps. For instance, looking at the expression 
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f'if'fj f'fj /) we observe that /'" must be a third derivative, since three arguments 
make it three-lmeav. This structure can be expressed in general by rooted trees, e.g., 



expressing /'/"(/./), 


expressing f"{f'f,f). 


Every node with n children denotes a nth derivative of /, which is applied as a 
multilinear map to further elementary differentials, according to the structure of 
the tree. We start reading off this structure by looking at the root. This defines a 
recursive procedure, if we observe the following: Having removed the root and its 
edges, a rooted tree /3 decomposes into rooted subtrees /3i,..., /3„ with strictly less 
nodes. The roots of the subtrees /3i,... ,/3„ are exactly the n children of /3’s root. 
This way a rooted tree j3 can be defined as the unordered list of its successors 

(2) (3 = [Pi, , Pn \) #/3 = 1 + #/?! + ... + #/3„. 

Here, we denote by #/3 the order of a rooted tree P, i.e., the number of its nodes. 
The root itself can be identified with the empty list, © = []. 

An application of this procedure shows for the examples above that /'(/"(/, /)) 
is expressed by [[©,©]] and /"(/'(/),/) is expressed by [[©],©]. The reader will 
observe the perfect matching of parentheses and commas. In general the relation 
between a rooted tree P = [Pi,..., /3„] and its corresponding elementary differential 
f^^\x) is recursively defined by 

(3) = f^-Hx) • [f^^^\x),..., . 

The dot of multiplication denotes the multilinear application of the derivative to 
the n given arguments. Due to the symmetry of the n-linear map the order 
of the subtrees Pi,... ,Pn does not matter, which means, that depends in a 
well-defined way on P as an unordered list only. 

From 0 = [] we deduce = /. Analogously, each of the recursive definitions 
in the following will have a well-defined meaning if applied to the single root 0 = [ ], 
mostly by using the reasonable convention that empty products evaluate to one and 
empty sums to zero—a convention that is also observed by most computer algebra 
systems. 


4. A SIMPLE PROOF OF BUTCHER’S THEOREM 

We are now in a position to calculate and denote the Taylor expansion of the 
continuous flow <!>’' in a clear and compact fashion. 

Lemma 1. Given f G the flow ^'^x expands to 

t -#/3 

(4) <^'^x = x+ ^ -—apf^^\x)+0{TP+'^). 

#I3<P 

The coefficients p\ and ap belonging to a rooted tree P = [Pi,..., /?„] are recursively 
defined by 

50 

«/3 = 


( 5 ) 


/3! = (#/3)/3i!-...-/3„!, 
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By 5f) we denote the number of different ordered tuples {/3i ,..., /3„) which correspond 
to the same unordered list (3 = [/3i,..., /3„]. 


Proof. The assertion is obviously true for p = 0. We proceed by induction on p. 
Using the assertion for p, the multivariate Taylor formula and the multilinearity of 
the derivatives we obtain 


#/3 


f{<^^x) = / I a; + + 0{tP+^) 


#0<P 


#/3i 


i: .^ 


n—0 


i#/3i<P 


#/3n<P 


r#/3n 


I + 0 {tP+^) 


■P+l'; 


r#/3l + ---+#/3n 


n=0 #/3i + ...+#/3„<p 

j(n).(/(/3,)^... j(/3n))+o(^P+l) 

= E E -- ^ ap, ■ ■ ■ ■ ■ f^^^ + 0{tP+ ) 


■i=0/3=[/3i....,/3„] 

#/3<P+l 


=ap 


#/3<P+l 

Plugging this into the integral form of the initial value problem we obtain 

l<T _ 

^'^x = x+ / f{^'^x)da = x+ 'Yj P 

#/3<p+l 


which proves the assertion for p + 1. 


□ 


A likewise clear and compact expression can be calculated for the Taylor expan¬ 
sion of the discrete flow. 

Lemma 2. Given f £ (7^(0, R"^) the discrete flow expands to 

( 6 ) ^^x = x+ Y f'^l^\x)+OiTP+^). 

#/3<p 

The vector £ R®, (3 = [fi,... ,(3A\, is recursively defined by 

(7) = (a-, i = l,...,s. 

Proof. Because of the definition (|^) of the discrete flow we have to prove that the 
stages ki expand to 

h = Y +0{tP). 

*P<P 

This is obviously the case for p = 0. We proceed by induction on p. Using the 
assertion for p, the definition of the stages fci, the multivariate Taylor formula and 
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the multilinearity of the derivatives we obtain 

ki = f (x + T I ^ + O(r^) 


K#P<P 


n—0 


“ ”■ \#/3i<P 

J2 (^A ■ A^^-^y +0(tP+1) 


#/3n<P 


“Hi n 


n=0 ”■ #/3i + ...+#/3„<p 


^#/3i+...+#/3„ .af3^-...-af3^yA- ■ 


... • (^ • ..., + 0(rP+i) 

= E E ^ «/3. ■ ■ • ■ • «/3„ ■ 4""^ + 0(r^>+i) 


n=0 /3=[/3i,...,/3„] 
#/3<p+l 


= Q;/3 


#/3<p+l 

which proves the assertion for p + 1. 


□ 


Comparing the coefficients of the elementary differentials in the expansions of 
both the phase flow and the discrete flow, we immediately obtain Butcher’s theo¬ 
rem 1^. 

Theorem 3 (Butcher 1963). A Runge-Kutta-method {b,A) is of order p £ N for 
all f £ CP{n,W^), if the order conditions 

(8) = 1//3! 

are satisfied for all rooted trees f3 of order fffd < p. 


5. Generating order conditions with Mathematica 

The recursive constructions underlying the proof of Butcher’s theorem can easily 
be realized using modern computer algebra systems like Mathematical We assume 
that the reader is familiar with this particular package. 

We begin by defining the recursive data-structure of a rooted tree as an unordered 
list—together with two simple routines for input and output: 

Attributes[Tree]={Orderless,bistable}; 

ToTree[f_]:=Tree@@ToTree/@f; ToTree[f_Symbol]:=Tree[]; 

Format [/3_Tree] ; =StringReplace [ToString [List@(§/3/. {}->"0"] , 

[" , "}"->"]"}] 

Now, here is an example for the input of the tree representing the elementary 
differential /"(/"(/, /'(/)), /) = /"(/, /"(/'(/), /)): 


Mathematica notebook with all the programs and examples of this e-print is included with 
the source files: http://arxiv.org/e-print/math.NA/0211049.tar.gz . 
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^i=ToTree [f " [f" [f , f' [f ] ] , f ] ] ; /32=ToTree [f" [f , f" [f' [f ] , f ] ] ] ; 

If [/3 i==/32./3i. .] 

[©,[©,[©]]] 

The definition (^) of the order #/3 is simply expressed by the recursive procedure: 
TreeOrder [/3_] : =l+Plus@@TreeQrder/0/3 

As an example, we take the elementary differential f'if'if'if), f'if), f), /): 

(3 = ToTree[f"[f'"[f'[f] ,f] ,f]] ; TreeOrder [/3] 

8 

The definition i) of /?! is analogously expressed by: 

TreeFactorial [/3_] : =TreeOrder [,9] Times@@TreeFactorial/@/3 
The above used tree [3 of order 8 gives: TreeFactorial [,9] 

192 

For the sake of completeness we also express the recursive definition (||) of ap using 
Mathematica: 


TreeAlpha[,9_] : =Length [Permutations [/3] ] /Length [/3] ! Times@@TreeAlpha/@/3 


An application to the above example yields: TreeAlpha[/3] 

1 

2 

We are now ready to map the recursive definition (^) of and of the order 
condition b'^ = 1//3! to Mathematica: 

TreeA[,9_,n_ : 1] : = (vars={i, j ,k,l,m,p,q,r,u,v,w}; 

Times@@(Sum[a„ars[[n]],vars[[n+i]]TreeA [#,n+l] //Evaluate, 

{vars [ [n+1] ] , s}//Evaluate] &/@,9) ) 

TreeOrderCondition [/3_] : =Sum[biTreeA [/3] //Evaluate, {i , s}] 
==l/TreeFactorial [,9] 

For convenience, the coordinate index of the vector A^^'^ can be chosen from the 
list {i,j,k,l,m,p,q,r,u,v,w} and is passed by number as the second argument to 
TreeA. The order condition belonging to the above example is obtained as follows: 

TreeOrderCondition [/3] 

S S 

'^j,k y]] O'k,! 

k=l 1=1 

Even the typesetting of this formula was done completely automatically, using 
Mathematica's ability to generate T^^^X-sources. 

To generate all the order conditions for a given order p, we need a device that 
constructs the set of all trees (3 with ^j3 < p. There are, in principle, two different 
recursive approaches: 




1=1 


1=1 


k=l 
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• root-oriented: generate all trees /3 of order ^/3 = phy first, listing all integer 
partitions p — 1 = pi + ... + pn, n = 1,... ,p — 1, and next, setting f3 = 
[/3i,..., /3„] for all trees /3i,..., /3„ of order #/3i = pi < p,... ^:^f3 = p^^ < p. 
These trees have already been generated by the recursion. 

• leaf-oriented: Add a leaf to each node of the trees /3 = [/3i,..., /3„] of order 
ff fd = p — 1, increasing thereby the order exactly by one. This can be done 
recursively by adding a leaf to every node of the subtrees /3i,..., 

The root-oriented approach was chosen by Sofroniou |P in his Mathematica package 
Butcher.m. It requires an efficient integer partition package and the handling of 
cartesian products. The leaf-oriented approach is as least as efficient as the other 
one, but much easier to code: 

AddLeave [/3_] : =AddLeave [/3] =Fold [Union, {Prepend[/3,Tree [] ] } , 

ReplacePart [/3, AddLeave [/3 [ [#] ] ] , #] &/ORange [Length [/3] ] ] 

Trees[order_] :=NestList[Union@@AddLeave/@#&,{Tree[]},order-1] 

Given an order p this procedure generates a list of the sets of trees for each order 
q < p, e.g.. Trees [4] 

{{©}, {[ 0 ]}. {[[©]], [©, ©]}, {[[[©]]], [[©, ©]], [©, [©]], [©, ©, ©]}} 

For instance, the number of trees for each order p < 10 is given by the entries of 
the following list: Length/QTrees [10] 

{1,1,2,4,9,20,48,115,286,719} 

The number of order conditions for p = 10 can thus be obtained by: PlusO®"/, 

1205 


Finally, for concrete calculations one has to specify the number s of stages. The 
following procedure then generates the specific set of equations for explicit Runge- 
Kutta methods: 


explicit={ai_.j_: >0/; i<=j , Ci->0} ; 

OrderConditions [order_, stages.] : = ( 

autonomous=Table[Sum [aij,{j,stages}]==ci, {i, stages}] ; 

{(TreeOrderCondition/QUnionOSTrees[order])/.s->stages/. 

ToRules[And@@autonomous].autonomous}/.explicit) 

This way, we can automatically generate and typeset the order conditions for the 
classical explicit 4-stage Runge-Kutta methods of order 4: 

First [OrderConditions[4,4] ] 

{ 6 i -F 62 + &3 + ^4 == 1; l>2 C2 + bsCs + b4 C 4 == —, 


bs C2 03,2 + ^4 (C2 04,2 + C3 04^3)--, 64 C2 03,2 04^3 == —, 


^3 ci 03^2 + ^4 (C2 04^2 + C3 04^3)-—, 62 C2 -F 63 C3 -F 64 C4--, 

bs C2 C 3 03,2 + bi C4 (C2 04,2 + C3 04,3) == &2 -F 63 C3 -F 64 ^4 == ^ } 
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Remark 4. Even in the more recent literature one can find examples like Q, 
where order conditions for Runge-Kutta methods are generated by using a computer 
algebra system to calculate the Taylor expansions of the flow and the discrete flow 
directly. This approach is typically bound to sco/ar non-autonomous equations, i.e., 
d = 1. Besides being inefficient for higher orders, it is well-known that for p > 5 
additional order conditions for general systems make an appearance, which do not 
show up in the scalar case. 

6. Examples of usage 


The following simple procedure tempts to solve the order conditions for a given 
order p and stage number s by using brute force, i.e.. Mathematicals Solve-com- 
mand. To simplify the task, the user is allowed to supply a set pre of a priori 
chosen additional equations and assignments that he thinks to be helpful. 

RungeKuttaMethod [p_, s_,pre_] : = ( 

rkTemplate={Table[aij,s}],Table[bi, , 

Table [ci,{j,1},{i,s}]}/.explicit; 
conditions=pre UUnionOOOrderConditions[p,s] ; 
solveVars=Complement[Flatten[rkTemplate],{0}] ; 
sol=Solve[conditions,solveVars]; 

Thread[{A,b,c}==MatrixForm/(§rkTemplate/.#]&/@sol) 


Since Mathematicats Solve-command uses a Grdbner-basis approach for solving 
systems of polynomial equations, we can show this way that the classical explicit 
4-stage Runge-Kutta method of order p = 4 is uniquely given by the additional 
constraints 62 = ^3 and C2 = C3: RungeKuttaMethod[4,4,{b2==b3,C2==C3}] 



(0 

0 

0 



1 

2 

0 

0 

0 

0 

1 

2 

0 

0 



0 

1 

0 / 



The next example is more demanding. In his book p. 199], Butcher describes 
an algorithm for the construction of explicit 6 -stage methods of order p = 5. The 
choices Cq = 1 and 62 = 0 together with the free parameters 02 , 03 , 04,05 and 043 
yield a unique method. Butcher provides a two-parameter example by choosing 
02 = M, 03 = 1 / 4,04 = 1 / 2,05 = 3 / 4,043 = V. By just passing this additional 
information to Mathematica's Solve-command we obtain the following solution: 


pre = {c2==u,C3==l/4,C4==l/2,C5==3/4,C6==l,b2==0,a4,3==v}; 
First[RungeKuttaMethod[5,6,pre]] 


/ 0 

u 


0 0 0 0 
0 0 0 0 


0 \ 

0 


{ 


A== 


— 1+8 u 
32 u 

— 1+4 r+2 v — ^uv 
8 u 


1 

32 u 
l-2v 
8 u 


0 000 

V 0 0 0 


3 (1—3 u—v-\-A uv) 3 ( —l+Ti) 

16 u 16 u 

— 7+22 r+6 v — 24uv 7—6 v 


14 u 


14 u 


-I (-1 + 0) ^ 0 0 

12 V 8 A 

7 7 7 ^/ 


b 


(J- n 1^ -2- IS 1^1 

V90 ^ 45 15 45 90 / ’ 


0 
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This result shows that the coefficients 051 and 052 of Butcher’s solution |^, p. 199] 
are in error, a fact that was already observed by Sofroniou ^ using the Mathematica 
package Butcher.m. 
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